# -*- coding: utf-8 -*-
"""
    @project: pythonProject
    @Author：HanYonghua
    @file： farrow_62.5_61.44_2.py
    @date：2025/7/8 20:05
    @blogs: https://www.ncatest.com.cn
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# 系统参数
fs_in = 62.5e6       # 输入采样率
fs_out = 61.44e6     # 输出采样率
ratio = fs_out / fs_in
impedance = 50       # 系统阻抗（欧姆）
ref_power_mw = 1.0   # dBm参考功率（1毫瓦）


# Farrow系数矩阵（保持不变）
cMatrix = np.array([
    [-8.57738278e-3, 7.82989032e-1, 7.19303539e+0, 6.90955718e+0, -2.62377450e+0, -6.85327127e-1, 1.44681608e+0, -8.79147907e-1, 7.82633997e-2, 1.91318985e-1, -1.88573400e-1, 6.91790782e-2, 3.07723786e-3, -6.74800912e-3],
    [2.32448021e-1, 2.52624309e+0, 7.67543936e+0, -8.83951796e+0, -5.49838636e+0, 6.07298348e+0, -2.16053205e+0, -7.59142947e-1, 1.41269409e+0, -8.17735712e-1, 1.98119464e-1, 9.15904145e-2, -9.18092030e-2, 2.74136108e-2],
    [-1.14183319e+0, 6.86126458e+0, -6.86015957e+0, -6.35135894e+0, 1.10745051e+1, -3.34847578e+0, -2.22405694e+0, 3.14374725e+0, -1.68249886e+0, 2.54083065e-1, 3.22275037e-1, -3.04794927e-1, 1.29393976e-1, -3.32026332e-2],
    [1.67363115e+0, -2.93090391e+0, -1.13549165e+0, 5.65274939e+0, -3.34847578e+0, -6.20715544e-1, 2.06619782e+0, -1.42159644e+0, 3.75075865e-1, 1.88433333e-1, -2.64135123e-1, 1.47117661e-1, -4.71871047e-2, 1.24921920e-2]
]) / 12.28


# Farrow重采样函数（保持不变）
def farrow_resample(signal, cMatrix, n_samples_out):
    n_samples_in = len(signal)
    order = cMatrix.shape[0] - 1
    n_taps = cMatrix.shape[1]
    
    input_index = np.arange(n_samples_out) / n_samples_out * n_samples_in
    input_index_integer = np.floor(input_index).astype(int)
    input_index_fractional = input_index - input_index_integer
    
    output_signal = np.zeros(n_samples_out)
    
    for ix_order in range(order + 1):
        x = input_index_fractional ** ix_order
        for ix_tap in range(n_taps):
            c = cMatrix[ix_order, ix_tap]
            data_ix = np.mod(input_index_integer - ix_tap, n_samples_in)
            delayed = signal[data_ix]
            output_signal += c * delayed * x
    
    return output_signal


# 修正后的dBm转换函数（关键修改）
def signal_power_to_dbm(signal_data, fs, nperseg=4096, impedance=50, window='blackman'):
    """
    计算信号的功率谱并转换为dBm单位
    
    参数:
    signal_data: 输入信号
    fs: 采样率
    nperseg: 每段FFT的长度
    impedance: 系统阻抗（欧姆）
    window: 窗函数类型
    
    返回:
    f: 频率轴
    psd_dbm: 以dBm为单位的功率谱
    """
    # 使用指定窗函数减少频谱泄漏
    f, psd_v2hz = signal.welch(signal_data, fs, nperseg=nperseg, window=window, return_onesided=True)
    
    # 转换为功率谱密度（W/Hz）
    psd_w_per_hz = psd_v2hz / impedance
    
    # 转换为毫瓦/Hz（mW/Hz）
    psd_mw_per_hz = psd_w_per_hz * 1000
    
    # 计算每频率bin的功率（mW）
    bin_width = fs / nperseg  # 频率分辨率
    psd_mw = psd_mw_per_hz * bin_width
    
    # 计算dBm：10×log10(功率(mW))
    psd_dbm = 10 * np.log10(psd_mw + 1e-15)  # 加小值避免log(0)
    
    return f, psd_dbm


# 抛物线插值法估计真实峰值
def parabolic_interpolation(spectrum, peak_idx):
    """
    使用抛物线插值法估计离散频谱的真实峰值位置和幅度
    
    参数:
    spectrum: 频谱数据
    peak_idx: 离散峰值的索引
    
    返回:
    x: 插值后的峰值位置（浮点数索引）
    val: 插值后的峰值幅度
    """
    if peak_idx > 0 and peak_idx < len(spectrum) - 1:
        y = spectrum[peak_idx-1:peak_idx+2]
        x = peak_idx + (y[2] - y[0]) / (2 * (2 * y[1] - y[2] - y[0]))
        val = y[1] - 0.25 * (y[0] - y[2]) * (x - peak_idx)
        return x, val
    return peak_idx, spectrum[peak_idx]


# 生成已知功率的测试信号
def generate_calibrated_signal(amplitude_v, freq_hz, fs, n_samples):
    """生成已知幅度的正弦波，用于验证dBm计算的准确性"""
    t = np.arange(n_samples) / fs
    signal = amplitude_v * np.sin(2 * np.pi * freq_hz * t)
    
    # 计算理论功率（dBm）
    rms_v = amplitude_v / np.sqrt(2)
    power_w = (rms_v ** 2) / impedance
    power_dbm = 10 * np.log10(power_w * 1000)
    
    print(f"测试信号参数：幅度={amplitude_v}V，频率={freq_hz/1e6}MHz")
    print(f"理论总功率：{power_dbm:.2f}dBm（基于{impedance}Ω阻抗）")
    
    return signal, power_dbm


# 主测试函数
def test_dbm_spectrum():
    # 生成已知功率的测试信号 - 增加样本数提高精度
    n_in = 16384  # 样本数增加到16384
    test_freq = 10e6  # 10MHz信号
    nperseg_value = 4096  # 增加FFT点数到4096，提高频率分辨率
    window_type = 'blackman'  # 使用Blackman窗减少频谱泄漏
    
    # 生成测试信号
    input_signal, input_power_dbm = generate_calibrated_signal(
        amplitude_v=1.0,  # 保持原始振幅
        freq_hz=test_freq,
        fs=fs_in,
        n_samples=n_in
    )
    
    # 重采样
    n_out = int(n_in * ratio)
    output_signal = farrow_resample(input_signal, cMatrix, n_out)
    
    # 计算频谱（dBm）
    f_in, dbm_in = signal_power_to_dbm(input_signal, fs_in, nperseg=nperseg_value, window=window_type)
    f_out, dbm_out = signal_power_to_dbm(output_signal, fs_out, nperseg=nperseg_value, window=window_type)
    
    # 计算理论峰值
    bin_width = fs_in / nperseg_value
    theoretical_peak = input_power_dbm - 10 * np.log10(bin_width) - 3.01  # -3.01dB来自单边谱修正
    
    print(f"频率分辨率：{bin_width:.2f}Hz（{nperseg_value}点FFT，{window_type}窗）")
    print(f"理论频谱峰值：{theoretical_peak:.2f}dBm")
    
    # 找出实际峰值位置并应用抛物线插值
    peak_idx_in = np.argmax(dbm_in)
    x_in, val_in = parabolic_interpolation(dbm_in, peak_idx_in)
    peak_freq_in = f_in[int(x_in)] / 1e6  # 转换为MHz
    peak_value_in = val_in
    
    peak_idx_out = np.argmax(dbm_out)
    x_out, val_out = parabolic_interpolation(dbm_out, peak_idx_out)
    peak_freq_out = f_out[int(x_out)] / 1e6  # 转换为MHz
    peak_value_out = val_out
    
    print(f"插值后频谱峰值：输入信号 = {peak_value_in:.2f}dBm @ {peak_freq_in:.6f}MHz")
    print(f"插值后频谱峰值：输出信号 = {peak_value_out:.2f}dBm @ {peak_freq_out:.6f}MHz")
    print(f"峰值误差：输入信号 = {peak_value_in - theoretical_peak:.2f}dB")
    print(f"峰值误差：输出信号 = {peak_value_out - theoretical_peak:.2f}dB")
    
    # 自动计算Y轴范围
    max_peak = max(np.max(dbm_in), np.max(dbm_out))
    y_min = min(-80, max_peak - 60)  # 确保至少有60dB的动态范围
    y_max = max_peak + 10  # 留出10dB的余量
    
    # 绘制频谱
    plt.figure(figsize=(14, 10))
    
    # 原始信号频谱
    plt.subplot(2, 1, 1)
    plt.plot(f_in / 1e6, dbm_in, label='原始信号', linewidth=1.5)
    plt.axhline(y=theoretical_peak, color='r', linestyle='--', label=f'理论峰值: {theoretical_peak:.2f}dBm')
    
    # 标记实际峰值
    plt.plot(peak_freq_in, peak_value_in, 'ro', markersize=8)
    plt.annotate(f'插值峰值: {peak_value_in:.2f}dBm\n@{peak_freq_in:.6f}MHz',
                xy=(peak_freq_in, peak_value_in),
                xytext=(peak_freq_in + 1, peak_value_in - 10),
                arrowprops=dict(facecolor='black', shrink=0.05, width=1.5, headwidth=8),
                fontsize=9)
    
    plt.title(f'原始信号频谱（{fs_in/1e6}MHz采样，{nperseg_value}点FFT，{window_type}窗）')
    plt.xlabel('频率（MHz）')
    plt.ylabel('功率（dBm）')
    plt.xlim(0, fs_in/2/1e6)
    plt.ylim(y_min, y_max)
    plt.grid(True)
    plt.legend()
    
    # 重采样后信号频谱
    plt.subplot(2, 1, 2)
    plt.plot(f_out / 1e6, dbm_out, label='重采样后信号', color='orange', linewidth=1.5)
    plt.axhline(y=theoretical_peak, color='r', linestyle='--', label=f'理论峰值: {theoretical_peak:.2f}dBm')
    
    # 标记实际峰值
    plt.plot(peak_freq_out, peak_value_out, 'ro', markersize=8)
    plt.annotate(f'插值峰值: {peak_value_out:.2f}dBm\n@{peak_freq_out:.6f}MHz',
                xy=(peak_freq_out, peak_value_out),
                xytext=(peak_freq_out + 1, peak_value_out - 10),
                arrowprops=dict(facecolor='black', shrink=0.05, width=1.5, headwidth=8),
                fontsize=9)
    
    plt.title(f'重采样后信号频谱（{fs_out/1e6}MHz采样，{nperseg_value}点FFT，{window_type}窗）')
    plt.xlabel('频率（MHz）')
    plt.ylabel('功率（dBm）')
    plt.xlim(0, fs_out/2/1e6)
    plt.ylim(y_min, y_max)
    plt.grid(True)
    plt.legend()
    
    plt.tight_layout()
    plt.show()


if __name__ == "__main__":
    test_dbm_spectrum()    